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We present a new formulation of the Einstein equations that casts them in an explicitly first order, flux-conservative, hyperbolic 
form. We show that this now can be done for a wide class of time slicing conditions, including maximal slicing, making it 
potentially very useful for numerical relativity. This development permits the application to the Einstein equations of advanced 
numerical methods developed to solve the fluid dynamic equations, without overly restricting the time slicing, for the first time. 
The full set of characteristic fields and speeds is explicitly given. 
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Introduction. For decades the standard approach [1] 
ONto numerical relativity has been based on a direct appli- 
^Vation of the 3+1 formulation of the Einstein equations 
by Arnowitt, Deser, and Misner [2]. This important con- 
&©ribution laid the foundation for most numerical work in 
^jthe field. As convenient as this formulation is for numer- 
"^-nical relativity, the structure of the equations is extremely 
uncomplicated and most of the work consisted in develop- 
ing "ad hoc" numerical codes for every case considered. 
,_^This is in contrast with the situation of modern Com- 
£■ — putational Fluid Dynamics, which deals with first order, 
^Dflux conservative, hyperbolic (FOFCH) systems of equa- 
^j tions, for which many advanced numerical methods have 
been developed [3] based on the particular mathematical 
^^structure of the equations [4] . 

After many pioneering attempts [5-7] it was shown 
^fcliat the full set of 3D Einstein equations could be put 
I in the FOFCH form [8], similar to hydrodynamics. This 
kj^evelopment allowed the same advanced numerical meth- 
ods used in hydrodynamics, such as conservative schemes 
and modern shock capturing methods, to be applied to 
the Einstein equations for the first time. In a number 
of tests involving spherical black holes [9,10], ID gen- 
eral relativistic hydrodynamics [11] and 3D gravitational 
waves [12], this formulation showed some of its strengths 
over the standard approach. However, the price to pay 
for this original formulation was that it required the re- 
strictive assumption of harmonic time slicing. Although 
this slicing condition is singularity avoiding, making it 
potentially useful for studies of strongly gravitating sys- 
tems, it is just barely so [13]. For this reason metric 
and curvature components can grow without bound near 
the singularity and, without some sort of horizon bound- 
ary condition [10,14,15], this slicing is not well suited for 
black hole spacetimes. 

In this Letter, we will show that this FOFCH formal- 



ism can now be extended to a rather wide class of slicing 
conditions suitable for many different strong and weak 
field applications, including the time honored choice of 
maximal slicing. This new development opens the door 
to the use of very mature numerical methods in numerical 
relativity for many problems of interest, including black 
hole spacetimes. It should also facilitate mathematical 
studies of the Einstein Equations, as there is a vast liter- 
ature on systems of hyperbolic conservation laws of this 
form (see, e.g., Refs. [3,4]) that can now be applied to 
general relativity. 

Evolution equations. It is well known that Einstein's 
field equations can be decomposed into two sets: the evo- 
lution system and the constraints. The energy and mo- 
mentum constraints contain no second time derivatives. 
The remaining equations form the evolution system. This 
is just a matter of choice, because an evolution equation 
plus a constraint leads to another evolution equation with 
the same physical solutions (the ones obtained from ini- 
tial data which satisfy the constraints). 

The standard choice for the evolution system is to take 
the space components of the Ricci tensor, but one could 
choose instead the space components of the Einstein ten- 
sor or any other combination obtained by using the con- 
straints in a suitable way. The freedom to choose the 
coordinate gauge allows one to complete the evolution 
system in many different ways, and this can also lead to 
many different systems of equations, each one with its 
own structure. Nevertheless, we know that the physical 
solutions of all these systems (the ones obtained from 
initial data which satisfy the constraints) are equivalent, 
although some systems will be better suited for numerical 
studies. Below, we will use special choices of variables, 
and make use of the constraints in the evolution system 
and the gauge choice to derive a system with the FOFCH 
structure. 
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The standard 3+1 ADM evolution system is given by: 

(d t - £p)fij = -2a Kij (1) 

and the evolution equations: 

(d t - £p)Kij = -a i;j + a [Rff + trK K {j - 2K ij K {j 

-4 4) ] , (2) 

where indices are raised with the inverse matrix Y J of 
the space metric. This system is first order in time, but 
second order in space. To obtain a system which is also 
first order in space, we will introduce auxiliary variables 
which correspond to the space derivatives, 

A k = d k ln a , B£ = 1/2 3^ , D kij = 1/2 8 klij . (3) 

Note that B, D or the shift [3 are not tensor quantities 
on the constant time hypersurfaces. Nevertheless, we will 
formally raise and lower indices with the three-metric jij . 

One could then simply insert these quantities into the 
standard ADM equations to obtain a first order system. 
However, doing so blindly does not bring any particu- 
lar advantage. A careful choice of variables will trans- 
form the equations into the FOFCH form that is es- 
pecially suited to mathematical analysis and numerical 
treatment. In particular, the evolution system can be 
written as 

dt Kij + 8 r [-(] r Kij + a\ r j } = a Stj . (4) 

where the terms \ k - are given by 

A£- =£>£• + 1/2 (A 3 +W 3 -D 3 l) 

+ l/2S k (Ai+2Vi- Av) , (5) 

and we have noted 

V k = D k r r - D r rk . (6) 

Sij is a source term involving only the fields themselves 
and not their derivatives: 

% = ~ ~ 2Ki k K kj + trK K^ 
+ 2/ a (K ir B/ + KjrB{) 

+ \D kri D kr 3 + v\Xij - r,-fcrr/ r 

- (2D kr k - A r )(D ijr + D jir ) 

+ Ai(Vj - 1/2 D^) + Aj(Vi - 1/2 D lk k ) . 

With this first order formulation, one also needs to 
evolve the space derivatives. The simplest way of do- 
ing so is just to take the time derivative of Eq. (3) and 
interchange the order of space and time derivatives: 



d t A k + d r [-(] r A k +aQS r k ] 

= (2B k r -atrsS r k )A r , (7) 
dtBj + d r [-/3 r B l ; + <*()< 

= (2B k r -atr.sS r k )B; , (8) 
d t D ki j + d r [-l3 r D kij + a S r k (Kij - Sij)] 

= (2B k r -atrsS r k )D rlJ , (9) 

where we have noted 

(d t - /3 r d r )lna = -aQ , (10) 
(5 t -/3 r 5 r )/3 i = -2a 2 Q i , (11) 
Sij = (Bij + Bji)/a , (12) 

so choices of Q and Q 1 will determine the gauge condi- 
tions. 

The quantity V k introduced by Eq. (6) is very interest- 
ing. One can compute its time derivative from (9) but 
we will use the momentum constraint to transform that 
equation into 

d t V k + d r [-(] r V k +a( S r k -trsS r k )} 
= aG° k + aA r (K r k - tr K S r k ) 

+ (D k : - S' k D r f) K r s + (25/ -atrsS r k )V r 
-2(D^-SlD^ r )(aK r s -B r s ) . (13) 

We will consider V k as an independent quantity to be 
evolved with Eq. (13). In that way, Eq. (6) becomes an 
algebraic constraint between V k and the spatial metric 
derivatives D k ij (the algebraic form of the momentum 
constraint). This is crucial to ensure the hyperbolicity of 
the evolution system. 

The set of Eqs. (1), (4), (7), (8), (9), and (13), together 
with the gauge Eqs. (10) and (11), has the special form 
we seek. The entire nonlinear system of the Einstein 
evolution equations can be written in the form of a first 
order system of balance laws as 

d t u + d i F i (u) = S(u) (14) 

which is familiar from many branches of physics. It is 
essential to point out that first spatial derivatives of the 
fields occur only through the flux term, and the source 
terms do not contain any derivatives. Note that we have 
not yet specified the time derivatives (10), (11) of the 
lapse nor the shift. Previously this form required the use 
of the harmonic slicing condition. This is one of the key 
results of this paper. 

Characteristic Fields. We want to study under which 
conditions the evolution system is strictly hyperbolic. 
This point is crucial to apply advanced CFD numerical 
methods in which (the principal part of) the system is to 
be diagonalized by writing it in terms of the eigenfields 
which propagate along characteristic surfaces with their 
own characteristic speed. 
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When dealing with a first order system, one must first 
choose a fixed space direction to discuss hyperbolicity 
by considering only space derivatives along the selected 
direction. We will take the k coordinate axis so, in the 
following, k is never a dummy index. 

^From Eqs. (1), (10), and (11), it is clear that the co- 
ordinate normal lines will be characteristic (with speed 
—j3 k ). The corresponding characteristic fields are the 
lapse, the shift and metric components, and 

A k , , B k \ , D kHj (k' ^ k) . (15) 

Also, the light cones are characteristic surfaces (speeds 
— f3 k ± a\/j kk ). The corresponding eigenfields are 

l kk (K ik , - 8ik .) + S k s\, + \\ k , (k' ± k) . (16) 

It is easy to see that both tr K and D kr are alge- 
braically independent from the fields listed above and 
can not be then recovered from them. These quantities, 
together with the vector V k and the set Q, Q 1 , A k , B k 
cannot be properly computed until one specifies a gauge 
condition to evolve the right-hand-side terms Q, Q 1 in 
(10,11). In that sense, the set of characteristic fields is 
not complete. 

A new class of gauge conditions. We are interested in 
invariant slicing conditions. This means that the space- 
time slicing provided by our coordinate condition must 
be invariant under any transformation of the space coor- 
dinates of every slice. We must use then slicing scalars, 
like a, Q or tr K and their proper time derivatives (note 
that the shift "vector" does not behave as a slicing vec- 
tor; it is a vector under time independent transformations 
only). We also want to use an algebraic condition, be- 
cause of their simplicity and low computational cost. If 
we restrict ourselves to zero order scalars, we can play 
only with a and we get either a geodesic slicing or one 
of its generalizations. If we allow also first order scalars, 
we get both Q and tr K . The most general homogeneous 
algebraic condition is then 

Q- f(a)trK = , (17) 

where / is an arbitrary function. The geodesic slicing is 
then included as a subcase with / = 0. The maximal 
slicing condition (tr K = 0) is included also as a limiting 
case when / diverges. The / = 1 case corresponds to the 
harmonic slicing. Another interesting case is the "1+log" 
slicing [16,17], obtained when / = 1/a; it mimics max- 
imal slicing near a singularity, when the lapse collapses 
to zero. It has been very effective in evolving black hole 
spacetimes in ID [16] and 3D [17]. The term "1+log" 
arises from the expression of a in terms of y/j that one 
obtains when integrating Eq. (17) in the eulerian case 
(zero shift vector). Note however that the invariance of 
Eq. (17) ensures that one can apply it to obtain the same 
slicing even with a nonzero shift vector. 



With these choices we get characteristic cones associ- 
ated to the eigenfields 

ftrK 1 kk +2( S kk -tr S1 kk ) T y7l kT ^ r r (/ + 0) , 

(18) 

with characteristic speeds 

-/3 k ±ay/fi*, (19) 

and we will call these gauge speeds because of their ex- 
plicit dependence on the slicing condition. Gauge speed 
coincides with light speed only in the harmonic case 
(/ = 1) and this may be considered to be the distinc- 
tive quality of harmonic slicing. Gauge speed becomes 
infinite for a maximal slicing, as one would expect with 
the maximal slicing elliptic equation. We note that in 
the maximal slicing case both Q and A k are determined 
through the maximal slicing elliptic equation, and decou- 
ple from the other quantities. They are computed sepa- 
rately because Eqs. (17) and (18) provide now only one 
independent eigenvector (trK). The discussion in this 
section applies then to the remaining set of equations 
which, as we show below, can be fully diagonalized. 

The vector V k and the set Q l , B k l , D kr still can not be 
obtained until one specifies a shift vector condition. The 
simplest choice is a zero shift vector (Q l = 0, Bij = 0) 
so that 

Vi, A k -fD k r r (20) 

are the remaining eigenfields (zero speed). But one also 
has other choices of shift that still permit the system 
to be fully diagonalized [18]. The use of such shifts in 
numerical studies will be reported elsewhere. 

It is clear that negative values of / will lead to imag- 
inary gauge speeds. Moreover, the set of eigenfields is 
complete only if / ^ 0. This means that the evolution 
system will be strictly hyperbolic iff / > 0. Note also that 
gauges with / < 1 will have poor singularity avoiding 
behavior because gauge speed would be lower than light 
speed. Therefore, cases with / > 1 will look more ap- 
pealing for most Numerical Relativity applications, and 
many such choices have already been shown to work well 
for numerically evolved black hole spacetimes [16,17]. 

As a first test of these ideas we have applied this new 
formulation of the Einstein equations to numerical stud- 
ies of spherical black hole spacetimes and compared re- 
sults to those obtained with the standard ADM formula- 
tion of the equations [19] for a number of slicing condi- 
tions, including harmonic, "1+log", and maximal. In all 
cases the new numerical techniques allow the black hole 
to be evolved farther by orders of magnitude and with 
much higher accuracy. As a dramatic example of the 
power of this approach, in Fig. 1 we show the error in 
the evolution of the apparent horizon mass for the stan- 
dard approach with staggered leapfrog compared with 
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results from our code that uses the new formulation and 
advanced TVD methods [3]. Both results use identical 
resolution (200 grid zones) and maximal slicing. The 
ADM code crashes at t = 150M after errors of more 
than 100% develop. The new approach can be run in- 
definitely with errors of a few percent. This formulation 
with advanced numerical techniques will be particularly 
important in the 3D case with black holes [17], apparent 
horizon boundary conditions [14,10,15], and gravitational 
waves [12]. Based on the encouraging results in spherical 
symmetry, a 3D code using these techniques is well under 
development, and results will be reported elsewhere. 
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FIG. 1. Comparison between the standard ADM approach 
and the new FOFCH formulation. Evolution using the ADM 
approach cannot be continued beyond t = 150M, while the 
FOFCH continues indefinitely with a fraction of the error. 

Conclusion. We have presented a powerful new first 
order, flux conservative, hyperbolic formulation of the 
Einstein equations that can be used with a wide class of 
gauge conditions, including maximal slicing. There are 
many advantages in this new formulation. First and fore- 
most, it allows the use of advanced numerical methods 
that dramatically improve the accuracy and stability of 
numerical studies of the Einstein equations. Also, given 
that the system will be diagonalized in terms of eigen- 
fields propagating along characteristic surfaces, one can 
specifically design boundary conditions for each of the 
fields, depending on the sign and value of their charac- 
teristic speeds (in fact, this is the idea behind many of 
the numerical methods). This is particularly important 
for the propagation of gravitational waves, as the tech- 
niques automatically identify and evolve the eigenfields 
that are propagating ("radiative variables"), as opposed 
to metric functions. This should also be crucial in the de- 
velopment of apparent horizon boundary conditions, as 
proper understanding of the inner boundary causal struc- 
ture is required. The new class of algebraic slicings also 
gives insight into the singularity avoiding properties of 
general slicings and will allow strongly gravitating space- 
times. The analysis of the gauge speeds allows the use 
of proper causality conditions, and the study of potential 



problems when introducing new families of slicing gauges 
and shifts. The study of the structure of the remaining 
source terms will be important when the sources drive 
the evolution. Finally, this formulation allows the grav- 
itational field to be treated on the same footing, with 
the same methods, as the equations of relativistic hydro- 
dynamics. This combination should lead to a powerful 
approach to full GR hydrodynamics. As this formalism 
is already developed in the full 3D case, there is no diffi- 
culty in applying it there. 

Acknowledgements. We would like to thank Larry 
Smarr for helpful discussions. This work is supported 
by the DGICyT of Spain under project PB91-0335. J.M. 
acknowledges a PFPI Fellowship from MEC of Spain. 
We also acknowledge the support of NCSA and NSF 
grants PHY94-07882, PHY/ASC93-18152 (Arpa supple- 
mented), and INT94-14185. 



[1] J. W. York, in Sources of Gravitational Radiation, edited 

by L. Smarr (Cambridge, 1979) 
[2] R. Arnowitt, S. Deser, and C. W. Misner, in Gravitation: 

An Introduction to Current Research, edited by L. Witten 

(John Wiley, New York, 1962). 
[3] R. J. Leveque, Numerical Methods for Conservation Laws 

(Birkhauser Verlag, Basen, 1992). 
[4] P. D. Lax, Hyperbolic Systems of Conservation Laws and 

the Mathematical Theory of Shock Waves, CBMF-NSF 

11 (SIAM, Philadelpia, 1973). 
[5] A. E. Fischer, J.E. Marsden, Comm. Math. Phys. 28, 1 

(1972). 

[6] L. Smarr, Ann. N. Y. Acad. Sci. 302, 569 (1977). 
[7] C. Bona and J. Masso, Phys. Rev. D 40, 1022 (1989). 
[8] C. Bona and J. Masso, Phys. Rev. Lett. 68, 1097 (1992). 
[9] J. Masso, Ph.D. thesis, Univ. of Balearic Islands, 1992. 
[10] C. Bona, J. Masso, and J. Stela, Phys. Rev. D 51, 1639 
(1995). 

[11] C. Bona, J. Ibanez, J. Marti, and J. Masso, in Gravitation 
and General Relativity: rotating bodies and other topics, 
Vol. 423 of Lecture Notes in Physics, F. Chinea (Ed.), 
Springer- Verlag, 1993. 

[12] P. Anninos, J. Masso, E. Seidel, W.-M. Suen, M. Tobias, 
in preparation. 

[13] C. Bona and J. Masso, Phys. Rev. D 38, 2419 (1988). 

[14] E. Seidel and W.-M. Suen, Phys. Rev. Lett. 69, 1845 
(1992). 

[15] P. Anninos, G. Daues, J. Masso, E. Seidel, W.-M. Suen, 
Phys. Rev. D 51, 5562 (1995). 

[16] D. Bernstein, Ph.D. thesis, Dept. Of Physics, University 
of Illinois Urbana-Champaign, 1993. 

[17] P. Anninos, K. Camarda, J. Masso, E. Seidel, W.-M. 
Suen, J. Towns, to appear in Phys. Rev. D (1995). 

[18] C. Bona, J. Stela, J. Masso, and E. Seidel, in General 
Relativity, R. RufRni and M. Keiser, eds., World Scien- 
tific, London, 1995. 

[19] D. H. Bernstein, D. W. Hobill and L. L. Smarr, in Fron- 
tiers in Numerical Relativity, C. Evans, S. Finn, and D. 
Hobill, eds., Cambridge, 1989. 



4 



